Spectrum of stochastic evolution operators: local matrix representation approach 
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A matrix representation of the evolution operator associ- 
ated with a nonlinear stochastic flow with additive noise is 
used to compute its spectrum. In the weak noise limit a per- 
turbative expansion for the spectrum is formulated in terms 
of local matrix representations of the evolution operator cen- 
tered on classical periodic orbits. The evaluation of perturba- 
tive corrections is easier to implement in this framework than 
in the standard Feynman diagram perturbation theory. The 
result are perturbative corrections to a stochastic analog of 
the Gutzwiller semiclassical spectral determinant computed 
to several orders beyond what has so far been attainable in 
stochastic and quantum-mechanical applications. 
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I. INTRODUCTION 

Any dynamical evolution that occurs in nature is af- 
fected by noise. In a neuronal system the noise might be 
comparable in magnitude to purported underlying deter- 
ministic dynamics; in celestial mechanics the degrees of 
freedom omitted from a particular set of equations may 
be accounted for by very weak noise. Our task here and 
in two preceding papers is to systematically account 
for the effects of noise on measurable properties such as 
dynamical averages |^ in classical chaotic dynamical sys- 
tems. 

The theory is also closely related to the semiclassical 
expansions based on Gutzwiller's formula for the trace 
in terms of classical periodic orbits [Q in that both are 
perturbative theories (in the noise strength or Ti) derived 
from saddlepoint expansions of a path integral containing 
a Cantor set of unstable stationary points (typically pe- 
riodic orbits) . The analogy with quantum mechanics and 
field theory has been made explicit in [0 where Feynman 
diagrams were used to find the lowest nontrivial noise 
corrections. Unfortunately like its quantum counterpart, 
the Feynman diagram method for stochastic dynamics 
quickly becomes unwieldy at higher orders; rather than 
applying it directly we turn the argument around and 



suggest that the more efficient recent approaches of g] 
and the present paper be applied to difficult perturba- 
tive problems of quantum mechanics and field theory. 

An elegant method, inspired by the classical pertur- 
bation theory of celestial mechanics, is that of smooth 
conjugations In this approach the neighborhood 

of each saddlepoint is flattened by an appropriate co- 
ordinate transformation, so the focus shifts from the 
original dynamics to the properties of the transforma- 
tions involved. An elementary example is the Ulam map 
f{x) — Ax{l ~ x) which is solved exactly by the trans- 
formation X — s\vl'{tt9 /2) leading to the piecewise lin- 
ear tent map f{9) — 1 — \1 — 26\. In general there is 
no such explicit solution, but the expressions obtained 
for perturbative corrections are much simpler than those 
found from the equivalent Feynman diagrams. Using 
these techniques, we were able to extend the stochas- 
tic perturbation theory to the fourth order in the noise 
strength. 

Fourth order should be sufficient for most realistic cal- 
culations, but does not provide enough information to 
determine the convergence properties of the expansion, 
or determine eigenvalues beyond the first few. In this 
paper we develop a third approach, based on construc- 
tion of an explicit matrix representation of the stochastic 
evolution operator. Numerical implementation requires 
a truncation to finite dimensional matrices, and is less el- 
egant than the smooth conjugation method, but for high 
expansion orders (here eighth, but higher orders seem 
quite feasible) and many eigenvalues it is currently un- 
surpassed. As with the previous formulations, it retains 
the periodic orbit structure, thus inheriting valuable in- 
formation about the dynamics. 

In the following sections we define the stochastic dy- 
namics and show how to obtain matrix representations, 
both globally and located on the periodic orbits, as an 
expansion in terms of the noise strength a. The matrix 
elements are obtained from derivatives of the dynamics 
computed around each periodic orbit. We give as a nu- 
merical example the quartic map considered in both pre- 
vious papers, although the approach is very general and 
is by no means restricted to one dimension, to maps, or 
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to Gaussian noise. We find that up to eighth order, the 
cumulants converge super-exponentially with the length 
of periodic orbit and the expansion is now shown to be 
accurate to larger values of a. 



II. THE STOCHASTIC EVOLUTION OPERATOR 
AND ITS SPECTRUM 

An individual trajectory in presence of additive noise 
is generated by iteration 



(1) 



where f{x) is a map, f„ a random variable with the nor- 
malized distribution p(^), and a parametrizes the noise 
strength. In what follows we shall assume that the map- 
ping f{x) is one-dimensional and expanding, and that 
the ^„ are uncorrelated. A density of trajectories 4>{x) 
evolves with time as 

= (-C o 0„) (y) = dxC{y,x)(l)n{x) (2) 



where C is the evolution operator 



C{y,x) = 6„{y - f{x)) 
SAx) ^ fsix- aOpiOd^ = (^) . (3) 



For a repeller the leading eigenvalue of the evolution 
operator yields a physically measurable property of the 
dynamical system, the escape rate from the repeller. In 
the case of deterministic flows, the periodic orbit theory 
yields explicit formulas for the spectrum of C as zeros 
of its spectral determinant Our goal here is to ex- 
plore the extent to which such methods are applicable to 
systems with noise and to quantum systems. In particu- 
lar, we are interested in exploring the dependence of the 
eigenvalues v{(j) of C on the noise strength parameter a. 

The eigenvalues are determined by the eigenvalue con- 
dition 



F(cr, via)) = det(l - C/iy{a)) = 



(4) 



where F{a, 1/z) = det(l — z£) is the spectral determin- 
ant of the evolution operator C. Computation of such 
determinants commences with evaluation of the traces of 
powers of the evolution operator 



tr 



zC 



- = Vc„z", C„=tr/:", (5) 

1 — zL 



n=l 



which are then used to compute the cumulants Qn 
Qn{C) in the cumulant expansion 

oc 

det(l-z£) = l-^Q„z", 

71=1 



(6) 



by means of the recursion formula 

Qn = — [Cn — Cn-lQl — • • • ClQn-l) 

n 

which follows from the relation 

/ °° n \ 
det(l - zC) = exp - ^ — tr£" . 



(7) 



(8) 



Our task is to compute the cumulants Qn- We start by 
introducing a matrix representation for C 



III. MATRIX REPRESENTATION OF 
EVOLUTION OPERATOR 

As the mapping f(x) is expanding by assumption, the 
evolution operator (Q) smoothes the initial distribution 
(j^ix). Hence it is natural to assume that the distribution 
4>n{x) is analytic, and represent it as a Taylor series, in- 
tuition being that the action of C will smooth out fine 
detail in initial distributions and the expansion of (t>n{x) 
will be dominated by the leading terms in the series. 

An analytic function g{x) has a Taylor series expansion 



9{x) = E 



m=0 



m! dy^ 



:9iy) 



Expanding C{y,x) in Taylor series in y enables us to 
rewrite traces of £" as 



dxdy C{y,x)C{x,y) 



m.m' 




tC{v, x) 



■u=0/ 



— r 'Cfit, y) 



u=0 



Following H.H. Rugh we now define the matrix 
(m,TO' = 0, 1,2, ...) 



Qyni' 



dx C{y, x) 



(9) 



y=o 



L is a matrix representation of C which maps the 
component of the density of trajectories 4>n{x) in to 
the y™ component of the density 0„+i(j/), with y = f{x). 
The desired traces can now be evaluated as traces of the 
matrix representation L, tr£" = trL" . As L is infinite 
dimensional, in actual computations we have to truncate 
it to a given finite order. The Feynman diagrammatic 
and the smooth conjugation methods developed in the 
preceding papers ||l|,p| require no such approximations. 
However, as we shall see below, for expanding flows the 
structure of L is such that its finite truncations give very 
accurate spectra. 

Our next task is to evaluate the matrix elements of L . 
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IV. WEAK NOISE EXPANSION OF THE 
EVOLUTION OPERATOR 

We have written the operator £ in (|^) in terms of the 
Dirac delta function, £{x',x) = J 6{x' — f{x) — a^)p{^)d£_, 
in order to emphasize that in the weak noise hmit the 
stochastic trajectories are concentrated along the deter- 
ministic trajectory x' = f{x). Hence it is natural to 
expand the delta function in a Taylor series in a 



£{x',x) = S{x'~f{x)) 

i-ar 



S^^\x' - fix)) / Cpim 



E 



where S^^-^y) — -^-^5{y) . This yields a representa- 
tion of the evolution operator centered along the deter- 
ministic trajectory, with the Perron-Frobenius operator 
5{x' — f{x)), and corrections given by derivatives of delta 
functions weighted by moments of the noise distribution 



C{x',x)^S{x'-fix)) + 

n=2 



i-ay 



■a„<5(")(x'-/(a;)). 

(10) 



In our numerical tests we find it convenient to assume 
that the noise is Gaussian, p{^) = /^/V27r. For the 
Gaussian noise all a„ moments are known, and the weak 
noise expansion of C is 



C(x' , x) 



1 



71=0 



i!2" 



/(^)) 



= 5(x'-/(.T)) + y<5(2)(a;' 



/(^)) 



+ 



^5(^)(a:'-/(x)) + ... . (11) 



The choice of Gaussian noise is not essential, as the meth- 
ods that we develop here apply equally well to any other 
peaked smooth noise distribution, as well as space de- 
pendent noise distributions p(x,^). In any case, as the 
neighborhood of any trajectory is nonlinearly distorted 
by the flow, the integrated noise is never Gaussian, but 
colored. 



V. LOCAL MATRIX REPRESENTATION OF 
EVOLUTION OPERATOR 

Traces of powers of the evolution operator are now 
also a power series in a, with contributions composed of 
^(m) j-jj-^^-j _ 2;^_|_j^) segments. The contribution is non- 
vanishing only if the sequence xi, 2:2, a;„, a;„+i = xi 



is a periodic orbit of the deterministic map f{x). Thus 
the series expansion of tr has support on all periodic 



points Xa 



of period n, f^^{xa) = Xa', the skeleton 



of periodic points of the deterministic problem also serves 
to describe the weakly stochastic flows. The contribution 
of the Xa neighborhood is best presented by introducing 
a coordinate system c/fa centered on the cycle points, to- 
gether with a notation for the map (^ and the operator 
(p[) centered on the a-th cycle point 



a = 1, 



., Hp 



/a(0) = f{Xa+4') 
£a(0a+l, 0a) = jO-iXa+l + (j)a+l,Xa + (j)a) ■ (12) 

The weak noise expansion ( p^ ) for the a-th segment op- 
erator is given by 



00 , , ^ 
n=0 



n 



a„<5(")(0'+^a+l-/a(0)). 



Repeating the steps that led to (Q) we construct the 
local matrix representation of Ca centered on the Xa — > 
Xa-i-i segment of the deterministic trajectory 



90'™' 



'^a(0', 4>) j- 

m! 



'=0 



V" (-<^)" 

^ — ^ 71! 

n— maa;(m — m',0) 



(13) 



Due to its simple dependence on the Dirac delta function, 
B can expressed in terms of derivatives of the inverse of 
/a(0): 



(Ba), 



Qn 



m.' 



del)' 

a" (f-\xa+l+<j^')~Xa) 



.'=0 



50'" m\\fi{fa\xa+^+4>'))\ 
(to -I- 1)! 



90 



In+l 



4>'=0 



'=0 



(14) 



where we introduced the shorthand notation J-'a{(t>') — 

fa^{Xa+l + 0') - Xa- 

If we expand Ta{4>') in ^ Taylor series, the constant 
term is zero, since fa^{xa+i) — Xa- So we can write: 



^a(0') 



00 (/) 
1=1 



(15) 



where 1/.F^'^ = /;;. 

The matrix elements can be calculated explicitly as a 
multinomial expansion p| 




°° j-n 



{n\ai, ...,a„)'a;°\..x°" 



(16) 
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where the second sum (^) goes over all non-negative 
integers such that: 

ai + 2a2 + ... + nan = n, oi + 02 + ... + a„ = to, (17) 
and the multinomial coefficient is: 

nl 

(nlai, 09, .... a„) — ^-r^ ttttt ; — -, — rr r. (18) 

^ ' ' ' ' ' (l!)°iai!(2!)'^2a2!...(n!)°"a„! ^ ^ 

We apply the formula ( p^ ) to Ta{4>') with power to + 1: 
(.F,(</>'))™+^-(TO+l)! ^ ^^(Z|ai,a2,...,a,)' 

• (.Fa(l))'^n-^a('^)"^..(.Fa('))'^'. (19) 

For the (n + 1) -th derivative of this expression evaluated 
at 0' = only the Z = n + 1 term is non- vanishing. 
The matrix elements vanish for n < m, so B is a lower 
triangular matrix: 

(Ba)„m = '^{n + l|ai, 02, a„+i)' 



where trLp = trL„j^L2 • • • Li is the contribution of the 
p cycle, and the power series in &■> follows from the ex- 
pansion (|l3| ) of Lq in terms of Bq. The subscript saddles is 
a reminder that this is a saddle-point approximation to 
tr£" (see ref. ||] for a discussion), valid as an asymptotic 
series in the limit of weak noise. 

As a simple check of the above formulas, consider the 
noiseless case, for which the (La)m'm = (Ba)m'm nia- 
trices are a representation of the deterministic Perron- 
Frobenius operator C\^^q. The are triangular with 
diagonal elements {^a)mm = ■ The trace of the C 

on a periodic orbit p is therefore 

trL,.trL,,L2...L,.|:^^^, 

and we recover the standard deterministic trace for- 
mula for the Perron-Frobenius operator 

oo ^ 

tr£" = ^nj,^<5„.„^,— — . (23) 

p r— 1 ' P' 



(20) 



The diagonal and the nearest off-diagonal matrix ele- 
ments can easily be worked out. Here we show the first 
four expressed in terms of the derivatives of the original 
map: 



a) m+1 ,m 
(Ba)m+2,m 



1 



_i(TO + 2)! r: 

2 



to! 



(21) 



24TO!|/^|/r V/, 

\fa\fa 



f/4 
J a 



48to! 



+ (to + 5)(to + 6) 



rl/ rill 

277r-4(- + 5)^ 

^ a J a 

J a 

fie 

J a 



where fa, fai ' ' ' refer to the derivatives of f{x) evaluated 
at the periodic point Xa- 

By assumption the map is expanding, |/^| > 1. Hence 
the diagonal terms drop off exponentially, as ^/\fa\"^~^^, 
the terms below the diagonal fall off even faster, and we 
are justified in truncating Bq, as truncating the matrix 
to a finite one introduces only exponentially small errors. 

In the local matrix approximation the traces of evolu- 
tion operators are approximated by 



tr£"| 



^ ^ Tip ^ ^ SrijUpr^^ Lp 



J=0 



(22) 



VI. PERTURBATIVE CORRECTIONS TO 
EIGENVALUES 

The eigenvalue condition (^ is an implicit equation 
for the eigenvalue v = i^{<j) of form F{a, i^(cr)) = 0. As 
the eigenvalue condition is satisfied for any a, all total 
derivatives of the eigenvalue condition with respect to a 
vanish, leading to 



d dvdF 



= 



da 
d^vdF 
d(j^ dv 
d^vdF 
da'^ dv 



da dv 

d^F 



OF 

^dv d^F 



dv^ 
^d^vdv d^F 
da^ da dv'^ 



,d^v d^F 
da'^ dadv 

dv d^F 



dv 
da 



da dadv 
dv\^ d^F 
da 
d^F 



d^F 
'da^ 



(24) 



dt 



dadv'^ 



d^F 



da da^dv da^ ' 

and so on. v{0) can be computed by cycle expansions for 
a deterministic, noiseless flow, a ^ then parametrizes a 
weak perturbation to the deterministic Perron-Frobenius 
operator C\^^q. The above formulas enable us to com- 
pute recursively, order by order in cr", the perturbative 
corrections to the eigenvalues of C 



v{a) = ^ Vmcr"^ , 



m— 



Id™ 

to! dcr"* 



, (25) 



(T = 



in terms of partial derivatives of the eigenvalue condition 
F{a,v{a)) 



4 



F, 



kl 



Qk+l 



(26) 



(T=0,1' = I^(0) 

In this notation the formulas (04) for take the form 



V2 = 



1 
1 

3!Fi 
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(^^02 + 2FiiJ.i + 2F2o 

(Foi + 3 F12 Fii V2 

+ 3 F21 + 6 F20 vx V2 + vf) 



(27) 



As shown in ref. F^i can be computed from exphcit 
cycle expansions. However, in numerical calculations we 
find it more expedient to proceede by first expressing the 
spectral determinant F in terms of the cumulants. The 
traces of L" evaluated by (|l^) yield a series in , and 
the coefficients Qnj in the cumulant expansion 



F = det(l - z/:) = 1 - 2^ 2^ Qnjz'^a^ (28) 

n=l j=0 

are then obtained recursively from the traces, as in (^: 

(29) 



Qnj — ~ 1 Cnj ^ ^ ^ ^ ^ Qk,j~lCn~k,l 
^ \ k=l l=Q 



This gives F = F{z = \/v ,a) and the partial derivatives 
Fki can be found. Substituted in ( |27| ) they yield the per- 
turbative corrections to the eigenvalues. The above cal- 
culations can be efficiently done by manipulating formal 
Taylor series. 



VII. NUMERICAL TESTS 

Here we continue the calculations of the eigenvalue cor- 
rections described in refs. where more details and 
discussion may be found. We test our perturbative ex- 
pansion on the repeller of the 1-dimensional map 



(30) 



This repeller is a clean example of an "Axiom A" expand- 
ing system of bounded nonlinearity and complete binary 
symbolic dynamics, for which the deterministic evolution 
operator eigenvalues converge super-exponentially with 
the cycle length 0. 

We start the numerical calculations by determining all 
prime cycles up to a given length. For each prime cycle 
p we compute the truncated evolution matrix Lp and its 
repetitions Lp*^ to the given order in cr, and evaluate the 
traces ( p2[ ) . For the map at hand we find that truncations 
of size [16 X 16] suffice to achive double precision accu- 
racy for most cycles. However, as the short orbits are less 



unstable, they require larger matrix truncations in order 
to attain the same precision, and we employ a [28 x 28] 
truncation for the 2-cycles, and a [34 x 34] truncation 
for the fixed points. With the coefficients in the traces 
expansion (^) evaluated numerically, the cumulants and 
the perturbative eigenvalue corrections follow from (|29| ) 
and (|2^). In case at hand, a good first approximation 
is obtained already at n = 2 level, using only 3 prime 
cycles, and n — Q (23 prime cycles in all) is in this ex- 
ample sufficient to exhaust the limits of double precision 
arithmetic. 




5 6 
n 

FIG. 1. The perturbative corrections ( p9| ) to the cumulants 
Inj plotted as a function of cycle length n (for perturbation 



orders j 
gence. 



0, 2, 4, 6, 8) all exhibit super-exponential conver- 



The size of the cumulants is indicated in fig. ^ and 
the perturbative corrections to the leading eigenvalue of 
the weak-noise evolution operator are given in table |. 
Encouragingly, the value of vq = 2076.47 . . . computed 
here is not wildly different to our previous numerical es- 
timate 1^ of 2700. Both the cumulants and the eigenvalue 
corrections exhibit a super-exponential convergence with 
the truncation cycle length n. The super-exponential 
convergence has been proven for the deterministic, 1^0 
part of the eigenvalue [0, but the proof has not been 
extended to the stochastic evolution operators. 

We have chosen to test the formalism on this simple 
example, as here we are in a fortunate situation that 
the escape rate for arbitrary noise strength a can be cal- 
culated numerically by other methods to a rather high 
accuracy. For example, one can discretize the stochastic 
kernel on a spatial lattice j^] and determine numerically 
the leading eigenvalue. 

The perturbative result in terms of periodic orbits and 
the weak noise corrections is compared to the eigenvalue 
computed by the numerical lattice discretizationin fig. ||, 
with the absolute difference between the numerical and 
the mth order perturbative results plotted. We see that 
the perturbative result z/(to, cr) = X^fc^o indeed 
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improves as more perturbative terms are added. 



10' 

10" 

10" 
v-v(m) 

10-^ 
10" 
10" 



10" 




0.01 



0.03 



0.05 



0.07 



0.09 



of techniques such as the Borel resummation, the asymp- 
totic expansions of general integrals of saddlepoint type, 
and asymptotics beyond all orders All of this is be- 
yond the scope of the present paper, and we defer a full 
discussion of asymptotics to a forthcoming paper pO||. 



IX. ACKNOWLEDGEMENTS 

G.V. and G.P. gratefully acknowledges the finan- 
cial support of the Hungarian Ministry of Education, 
FKFP 0159/1997, OMFB, OTKA T25866/F17166. G.V. 
thanks Bruno Eckhardt the cordial hospitalty at the De- 
partment of Physics of the Philipps-Universitat Mar- 
burg and the Humboldt Fundation for support. G.P. 
thanks the EU network "Pattern formation, noise and 
spatio-temporal chaos in complex systems", TMR con- 
tract ERBFMRXCT960085, for partial support. N.S. is 
supported by the Danish Research Academy Ph.D. fel- 
lowship. 



FIG. 2. The difference between the numerical and pertur- 
bative eigenvalue |i^(o-) — u{'m,a)\. The plateau at 10~^ is a 
numerical artifact due to the limited accuracy of the lattice 
discretization calculation. 



VIII. SUMMARY AND OUTLOOK 

In this paper we study evolution of a classical dynami- 
cal system with additive noise. In the limit of weak noise 
the traces of the corresponding evolution operator are ap- 
proximated by sums of local traces computed on periodic 
orbits. Here we present a new, computationally efficient 
technique for evaluation of these local traces based on 
a matrix representation of the evolution operator, and 
show that method is powerful enough to enable us to 
compute 2 more orders of perturbation theory. 

The local matrix representation can be interpreted as 
follows. Substituting (22|) into (||) we obtain 



det(l - z£)|_ = n det(l - ^"^Lp) ■ 



(31) 



In other words, in the saddle-point approximation the 
spectrum of the global evolution operator C is in this ap- 
proach pieced together from the local spectra computed 
cycle-by-cycle on neighborhoods of individual prime cy- 
cles with periodic boundary conditions. Vattay ||] was 
first to formulate the h corrections to the semi-classical 
Gutzwiller theory in terms of local spectra. Here we have 
shown that also the stochastic flows can be suspended on 
the skeleton of classical periodic orbits in this way. 

With so many orders of perturbation theory, we are 
now poised to address the issues raised by the asymp- 
totic series nature of perturbative expansions. We can 
now hope to resum the series to all orders, making use 
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n 




1^2 


1/4 






1 


0.308 


0.42 


2.2 


17.4 


168.0 


2 


0.37140 


1.422 


32.97 


1573.3 


112699.9 


3 


0.3711096 


1.43555 


36.326 


2072.9 


189029.0 


4 


0.371110995255 


1.435811262 


36.3583777 


2076.479 


189298.8 


5 


0.371110995234863 


1.43581124819737 


36.35837123374 


2076.4770492 


189298.12802 


6 


0.371110995234863 


1.43581124819749 


36.358371233836 


2076.47704933320 


189298.128042526 



TABLE I. Significant digits of tlie leading deterministic eigenvalue 1^0, and tlie a^,---,a^ perturbative coefficients (ps]), 
calculated from the cumulant exapansion of the spectral determinant, as a function of the cycle truncation length n. Note the 
super-exponential convergence of all coefficients. 
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